Kinetic Simulations of Magnetized Turbulence in Astrophysical Plasmas 
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This letter presents the first ab initio, fuUy electromagnetic, kinetic simulations of magnetized 
turbulence in a homogeneous, weakly coUisional plasma at the scale of the ion Larmor radius (ion 
gyroscale) . Magnetic and electric-field energy spectra show a break at the ion gyroscale; the spectral 
slopes are consistent with scaling predictions for critically balanced turbulence of Alfven waves above 
the ion gyroscale (spectral index —5/3) and of kinetic Alfven waves below the ion gyroscale (spectral 
indices of —7/3 for magnetic and —1/3 for electric fluctuations). This behavior is also qualitatively 
consistent with in situ measurements of turbulence in the solar wind. Our findings support the 
hypothesis that the frequencies of turbulent fluctuations in the solar wind remain well below the ion 
cyclotron frequency both above and below the ion gyroscale. 
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Introduction. A wide variety of astrophysical 
plasmas — e.g., the solar wind and corona, the in- 
terstellar and intracluster medium, accretion disks 
around compact objects — are magnetized and turbulent. 
The turbulence in these systems is damped at small 
scales where the plasma is only weakly collisional, so 
a kinetic description is required. It is often a good 
approximation to consider small fluctuations occurring 
on top of an equilibrium state with a uniform (or 
large-scale) dynamically strong mean magnetic field 
(the Kraichnan hypothesis 1]). The resulting (subsonic) 
magnetohydrodynamic (MHD) turbulence is believed to 
be a Kolmogorov-like cascade of spatially anisotropic 
Alfvenic fluctuations 0]. Such anisotropy is observed in 
laboratory plasmas 0] , the solar wind 0] , and numerical 
simulations Assuming a critical balance between 

the linear frequencies and nonlinear decorrelation rates 
0) 01 1 the anisotropy is scale-dependent with wave 
numbers parallel and perpendicular to the local mean 
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field related by fc|| cx kj^ . This implies that in most 
astrophysical plasmas, the frequencies of the Alfvenic 
fluctuations remain below the ion cyclotron frequency, 
to = k\\VA ^ J^i, even as the perpendicular wavelength 
reaches the ion gyroscale, k±pi 1. 

Such fluctuations are well described by gyrokinet- 
ics (GK), a rigorous low- frequency anisotropic limit 
of kinetic theory 0, H, [13]) which systematically 
averages out the cyclotron motion of particles about 
the magnetic field. In GK, the MHD fast wave and 
cyclotron resonances are ordered out, while finite 
Larmor radius effects and the collisionless Landau 
resonance are retained. GK enables numerical studies of 
kinetic turbulence with today's computational resources 
because the gyroaveraging eliminates fast time scales 



and reduces the dimensionality of phase space from six 
to five (three spatial dimensions plus the parallel and 
perpendicular particle velocities). GK has been used 
to study electrostatic turbulence in fusion plasmas for 
decades, but there have been few GK treatments of 
astrophysical plasma turbulence. GK is not applicable 
to large-scale, roughly isotropic fluctuations, such as 
are directly excited in the interstellar medium by 
supernovae. However, the fluctuations in magnetized 
plasma turbulence become smaller amplitude and more 
anisotropic at smaller scales. GK theory and simulations 
are thus appropriate, and hold considerable promise, 
for studies of microscopic phenomena such as turbulent 
heating and magnetic reconnection, and for interpreting 
observations of short-wavelength turbulent fluctuations. 
This Letter reports the first ab initio, fully electromag- 
netic, kinetic simulations of turbulence in a magnetized 
weakly collisional astrophysical plasma. 

The study of turbulence in weakly collisional plasmas 
benefits greatly from access to a unique laboratory — the 
near-Earth solar wind — in which spacecraft make in 
situ measurements of the properties of turbulence 
from the large (energy-containing) scales to the small, 
kinetic scales at which fluctuations are damped. The 
one-dimensional frequency spectrum of magnetic fluc- 
tuations typically shows a power-law behavior with a 
—5/3 slope at low frequencies [llj, a break at a few 
tenths of a Hz, and a steeper power-law at higher 
frequencies with a slope that varies between —2 and 
—4 [l§|. It is generally agreed that the —5/3 range 
is an MHD inertial range, while the break and the 
dissipation-range slope have been variously attributed 
to proton cyclotron damping Il3|, Landau damping of 
kinetic Alfven waves (KAW) or the dispersion of 
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whistler waves [15| . Recent simultaneous measurements 
of the magnetic- and electric-field fluctuations found an 
increase in the wave phase velocity above the spectral 
break [fij . a finding consistent with the conversion to a 
KAW cascade but inconsistent with cyclotron damping 
[lo| . The GK simulations presented below capture 
all of the spectral features described above and show 
magnetic- and electric-e nerg y spectra similar to those 
reported empirically in |f6l |. Our simulation results 
suggest that the turbulent fluctuation spectra observed 
in the solar wind are a consequence of the transition 
from an Alfven-wave to a KAW cascade. 

The Code. We have used AstroGK, a new GK code 
developed specifically to study astrophysical turbulence. 
AstroGK is essentially a slab version of the publicly avail- 
able code GS2, used to study plasma turbulence in fusion 
devices [l3|. We now give a brief overview of the code. 
A detailed description will appear elsewhere. 

The simulation domain is a periodic flux tube with a 
straight uniform mean magnetic field Bq and no equilib- 
rium gradients. The equilibrium distribution is taken to 
be Maxwellian for all particle species. The code solves 
the GK equation [sl, evolving the perturbed gyroaver- 
aged distribution function hs{x,y, z,es,^) of the guid- 
ing centers for each species s — ions (protons) and elec- 
trons with the correct mass ratio nii/me = 1836. Spa- 
tial dimensions (x, y) perpendicular to the mean field are 
treated pseudospectrally; a conservative finite-difference 
scheme is used in the parallel direction z. A gyroaveraged 
pitch-angle-scattering collision operator [9] is used. The 
pitch-angle derivatives are done using second-order fi- 
nite differences. The electromagnetic field is represented 
by the scalar potential ip, parallel vector potential 
and the parallel magnetic field perturbation 6B\\ . These 
are determined from the quasineutrality condition and 
Ampere's law 8], where the charge densities and cur- 
rents are calculated as velocity-space moments of the per- 
turbed distribution function. These velocity-space inte- 
grals (over particle energies Eg — and pitch angles 
^ = fy/w) are done with spectral accuracy, using high- 
order Gaussian-Legendre integration rules. The linear 
terms in the GK system, including the field equations, 
are advanced implicitly in time; for the nonlinear terms, 
an explicit, third-order Adams-Bashforth scheme is used. 

Linear Benchmarks. In an earlier paper [8,], we ver- 
ified that GS2 correctly describes linear kinetic physics 
in the parameter regimes relevant to astrophysical plas- 
mas. AstroGK has been checked to agree with GS2 ex- 
actly and also benchmarked against linear kinetic the- 
ory, as illustrated by Fig. [TJ for k±pi ^ 1, we have 
Alfven waves, u — ±k\\VA, and the damping is very 
small; for k±pi ^ 1, these become kinetic Alfven waves, 
Lu — ±k\\VAk^Pi/ \/ Pi -f 2/(1 -I- Te/Ti), so their phase ve- 
locity increases linearly with k± and they are also more 
strongly damped, in excellent agreement with linear the- 
ory [3, ll8|. Here va = Bq / y/AwrnjjTi is the Alfven speed. 



: ' 1 


"/(l<||Vj 






7/(k,vJ^ 


Theory -i 






■= AstroGK : 


". 1 


1 1 1 1 1 


/?,= 1 VT,= 1 1 
... 1 1 ' 


0.1 




1 10 



FIG. 1: Normalized frequencies u/k\^VA and damping rates 
'y/k\\VA vs. the normalized perpendicular wave number k±pi 
for a plasma with Pi = 1 and Ti/Te. — 1. AstroGK (open 
squares) correctly reproduces the analytic results from the 
linear coUisionless gyrokinetic dispersion relation (line) 0]. 

Ui the ion number density, T^. and Ti the ion and electron 
temperatures, and Pi — SnniTi/B^. 

Driving. The driving and dissipation scales in astro- 
physical turbulence are widely separated: e.g., in the slow 
solar wind, the ion gyroscale is pi ^ 10^ cm, while the 
effective driving scale is L ~ 10^"'^ cm Td]. Such scale sep- 
arations are, of course, not accessible numerically. In our 
simulations, the size of the domain is understood to be 
much smaller than the driving scale. We model the en- 
ergy influx from larger scales by adding to Ampere's law 
a parallel "antenna" current j^^^^. For each chosen driving 
wave vector , the antenna amplitude is calculated from 
a Langevin equation whose solutions arc Alfven waves 
with wave vector k^, frequency lo — ±fc^||W^ and a decor- 
relation rate comparable to lo. This method of driving is 
motivated by the theoretical expectation that the turbu- 
lence in the inertial range (at scales pi ^ fc^^ ^ i) is 
Alfvenic and critically balanced . 

Dissipation. The driving injects power into the sys- 
tem in the form of electromagnetic fluctuations. In 
steady state, this power must be dissipated into heat. By 
Boltzmann's H-theorem, no entropy increase and, there- 
fore, no heating is possible in a kinetic system without 
collisions. If the collision rate is smaller than the fluctua- 
tion frequencies, the perturbed distribution function de- 
velops small-scale structure in velocity space 0, 0] . This 
makes the velocity derivatives in the collision integral 
large so the collisions can act, a situation analogous to 
the emergence of small spatial scales in neutral fluids with 
small viscosity (Kolmogorov cascade). In GK turbulence, 
the cascades in position and velocity space are linked, so 
we may speak of a kinetic cascade in five-dimensional 
phase space 9]. CoUisionless Landau damping of the 
electromagnetic fiuctuations leads to particle heating in 
the sense that it transfers the electromagnetic fluctua- 
tion energy into fluctuations of the particle distribution 
function (the kinetic entropy cascade ^]), which are then 
converted into heat by collisions. 

A detailed analysis of the kinetic cascade will be pre- 
sented in a separate study, but the lesson is that ki- 
netic turbulence simulations need to include collisions 
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FIG. 2: (Color online) Magnetic (solid line) and electric 
(dashed line) energy spectra in the MHD regime {k±pi < 1). 
The box size is L±/2Tr = lOpi. Electron hypercoUisionality is 
dominant for k±pi > 1 (dotted line). 




and need to have sufficient velocity-space resolution for 
the correct relationship to be established between small- 
scale structures in velocity and position space. Accom- 
plishing this with a physical collision operator simultane- 
ously for ions and electrons is very difficult. To ease the 
resolution requirements, we employ a hypercoUisionality 
(analogous to hyperviscosity in fluid simulations). This 
takes the form of a pitch-angle-scattering operator with 
a wave- number-dependent collision rate i^ft,(A;_L/fcj_max)^, 
where fc^,„ax is the grid-scale wave number. This artifi- 
cially enhanced collision term terminates the cascade and 
produces positive-definite heating close to the grid scale, 
while allowing essentially coUisionless physics at larger 
scales. For the ions, the importance of the hypercoUi- 
sionality is marginal, while for the electrons, we needed 
a large value of Vh- As a result, electron heating (at the 
electron gyroscale pe) is not well modeled, but this is an 
acceptable sacrifice because our focus is on the turbulent 
cascade through the ion gyroscale at pe- 

Results. The physical parameters in GK simulations 
of plasma turbulence are Ti/Tg and j3i. Here both 
are set to 1, sensible characteristic values for the so- 
lar wind at 1 AU, and for the interstellar medium; a 
full parameter scan is clearly desirable in the future 
(e.g., (3i in the solar wind at 1 AU varies roughly be- 
tween 0.1 and 10). By varying the driving wave num- 
ber ka and the (hyper)collision rate, we may focus on 
various scale ranges. Here we present results obtained 
for the inertial range {k±pi <^ 1) and around the ion 
gyroscale {kj_pi ~ 1). In what follows, the normal- 
ized magnetic-energy spectrum is defined Eg^^k^) — 
{L,/L\)2'Kk^^J dz{\Aik^{z)\'^)/8T:n,T„ where k± is 
measured in units of p~^, and Lj_ are parallel and 
perpendicular box dimensions, and the angle brackets de- 
note angle averaging over a wavenumber shell centered at 
|kx| — k± and with the width equal to 2Tr/L±. The nor- 
malized electric-energy spectrum EEj_{k±) is defined in 
a similar way in terms of ipkj_ , with an extra factor of 
(c/va)^, where c is the speed of light. 

In the inertial range, k±pi ^ 1, the Reduced MHD 
equations are the rigorous limit of GK for Alfvenic fluc- 



FIG. 3: (Color online) Bold lines: normalized energy spectra 
for SB± (solid), J-By (dash-dotted), and -Ex (dashed). Thin 
lines; solution of the turbulent cascade model of [lO|. The 
resolution of this simulation is {N:,:, Ny, N^, N^:, N^, Ns) = 
(64,64,128,8,64,2), requiring ~ 0.5 x 10^ computational 
mesh points. The box size is L±/2tt — 2.5pi. Electron hyper- 
coUisionality is dominant for k±pi > 8 (dotted line). 

tuations [Oj]. Thus, kinetic turbulence in this regime must 
be consistent with the numerical results obtained in MHD 
simulations Q . Fig. [2] shows the normalized magnetic 
and electric energy spectra calculated gyrokinetically in 
this regime. As expected for critically balanced Alfvenic 
turbulence Q, these spectra are coincident and have a 
scaling consistent with k^^^^. This is the first demon- 
stration of an MHD turbulence spectrum in a kinetic 
simulation. While this is not a surprising result, it can 
be viewed as a fully nonlinear benchmark. 

Our main numerical experiment focuses on scales 
near k±pi ^ 1. This regime cannot be simulated by 
any fluid model. However, we know from theory that 
low-frequency Alfvenic turbulence is rigorously described 
by Reduced MHD equations for k±pi ^ 1 and by a 
similarly reduced version of the Electron MHD equa- 
tions for k±pi ^10. The latter system supports ki- 
netic Alfven waves (see Fig. [T]) . If one assumes a turbu- 
lent cascade of KAW-like fluctuations decorrelating on a 
timescale comparable to the linear KAW period (critical 
balance), a Kolmogorov-style scaling argument predicts 
that the magnetic-energy spectrum steepens from k^^^^ 

— 7/3 

to fcj^ , while the electric-energy spectrum flattens to 

A:^^^'^ 0, [13, [l^l- Thus, a spectral break is expected 
around k±pi ~ 1, corresponding to the transition be- 
tween Alfven-wave and KAW turbulence. Fig. [3] shows 
the energy spectra in our simulations near this transition. 
A spectral break is, indeed, observed (at k±pi ~ 2), as 
is the steepening (flattening) of the magnetic- (electric- 
)energy spectra. The spectra at wave numbers below 
and above the transition are consistent with the above 
predictions for critically balanced Alfven-wave and KAW 
cascades @,S,[i3,[ii- 

There is a striking similarity between the simulated 



4 



spectra shown in Fig. [3] and the magnetic- and electric- 
energy spectra in the solar wind reported in [l^ . The in- 
crease in phase velocity in the dissipation range {k±pi > 
1), shown by both measurement and simulation, is com- 
pelling evidence that the observed breaks in the spectra 
are caused by a transition to a KAW cascade, not by the 
onset of ion cyclotron damping 10]. 

The scaling predictions for KAW turbulence are made 
assuming negligible Landau damping. In our simulations, 
the damping is, indeed, small, so it is reasonable that the 
scaling predictions are well satisfied. However, this will 
not be true in all real astrophysical situations. We have 
argued that the spectra much steeper than k^^^ 
often observed in the solar wind IJ] can be due to non- 
negligible Landau damping. A simple way to estimate 
the effects of the damping on the energy spectra was 



proposed in [10| (see also [18|), where a spectral model 



of the turbulent cascade was constructed based on three 
assumptions: (i) spectrally local energy transfer, (ii) crit- 
ical balance, (iii) the applicability of the linear damping 
rates. Using this model, the energy spectrum EB^{k±) 
can be predicted in the entire simulation range, given 
one "Kolmogorov" constant, which quantifies the linear 
damping rate relative to the non-linear cascade rate. In 
Fig. [31 we show that this analytical model reproduces the 
entire shape of the numerical spectrum. The model works 
well without fine tuning, for a range of values of the con- 
stant; this is because the damping is small in this simula- 
tion and our model captures the transition from Alfvenic 
to KAW turbulence. The agreement between the ana- 
lytical model and the simulations is a non-trivial result: 
it suggests that the linear damping rate does not signif- 
icantly underestimate the rate at which the electromag- 
netic energy is dissipated in the nonlinear simulations. 
Future simulations will determine whether stronger lin- 
ear damping can account for the steeper spectra often 
observed in the solar wind. 

A further test of the conclusion that we are seeing a 
KAW cascade in the simulations is achieved by using the 
linear GK cigcnfunctions for KAWs to produce the en- 
ergy spectra for the electric fluctuations {E±) and for 
the fluctuations of the field strength (SB^^). These fit the 
spectra measured from our numerical data well (Fig. [3]). 

Conclusions. We have presented first-of-a-kind ki- 
netic simulations of turbulence in a weakly coUisional, 
magnetized plasma. The ion-gyroscale turbulent fluc- 
tuations simulated here represent the fate of a larger- 
scale MHD cascade. The qualitative agreement between 
our simulations and solar-wind measurements [l6l | sup- 
ports theoretical models in which the turbulent fluctua- 
tions in the solar wind have frequencies well below the 
ion cyclotron frequency even when the cascade reaches 
the (perpendicular) scale of the ion Larmor radius. The 
observed break in the magnetic-energy spectrum in the 
solar wind is inferred to correspond to a transition to 
kinetic- Alfven- wave turbulence, not to the onset of ion 



cyclotron damping. Although half a billion mesh points 
were used in the case of Fig. |3l the resolution in ve- 
locity space is still not fully sufficient to draw detailed 
conclusions about the turbulent heating. Nonetheless, 
the agreement between the simulations and an analyti- 
cal cascade model based on linear damping rates implies 
that the latter do not significantly underestimate the true 
damping in a turbulent coUisionless plasma. Future sim- 
ulations will probe a range of plasma parameters, includ- 
ing more heavily damped regimes, that will allow a more 
quantitative study of the role of coUisionless damping in 
turbulent plasmas. The first results reported in this Let- 
ter demonstrate that such kinetic simulations of plasma 
turbulence may be undertaken with some confidence, us- 
ing existing computational resources. 
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